#load pakages
if(!require(pacman)){install.packages("pacman")}
pacman::p_load(tidyverse, haven, ggpubr, rstatix, glue, estimatr, ggplot2, MetBrewer)

#load data
data <- read_dta("data/hgopy_pilot_public.dta")
#create the new price variable with treatments as factors
data <- data %>% mutate(data, 
                        pdx = ifelse(price_larc_cat==3,'Full','Discounted'),
                        pdn = ifelse(price_larc_cat==3,0,1), 
                        pd = factor(pdx),
                        pd = fct_reorder(pd,pdn),
                        csx = ifelse(view==0,'MiM',
                                     ifelse(view==1,'IDM',
                                            ifelse(view==2,'SDM','None'))),
                        cs = factor(csx),
                        cs = fct_reorder(cs,view))
#create the new SARCs variable
data <- data %>%
  mutate(
    adopted_sarc_alt1 = case_when(
      (adopted_sarc == 1) ~ 1,
      (adopted_other < 4) ~ 1,
      TRUE ~ 0
    )
  )
table(data$adopted_sarc)
table(data$adopted_sarc_alt1)
table(data$adopted_sarc_alt1,data$adopted_sarc)
table(data$adopted_sarc_alt1,data$adopted_other)
#create a LAM variable, take the first consultation only? 
data <- data %>%
  mutate(
    adopted_lam = ifelse(adopted_lam1==1,1,0),
    adopted_lam = ifelse(is.na(adopted_lam1),0,adopted_lam),
    adopted_lam = ifelse(adopted_larc==1,0,adopted_lam),
    adopted_lam = ifelse(adopted_sarc==1,0,adopted_lam),
  )
table(data$adopted_lam)
table(data$adopted_lam,data$adopted_lam1) 
table(data$adopted_lam,data$adopted_method) 
#now, a general sarc that also includes LAM
data <- data %>%
  mutate(
    adopted_sarc_alt2 = case_when(
      (adopted_sarc_alt1 == 1) ~1,
      (adopted_lam == 1) ~ 1,
      TRUE ~ 0
    )
  )
table(data$adopted_sarc_alt2)
table(data$adopted_sarc_alt2,data$adopted_sarc_alt1)
#then, adopted nothing at all... 
data <- data %>%
  mutate(
    nothing = case_when(
      (adopted_sarc_alt2 == 1) ~ 0,
      (adopted_larc == 1) ~ 0,
      TRUE ~ 1
    )
  )





###Plot - definition 1, broader definition SARC measure

###Plot prep 

# define outcome for convenience
data = data %>% mutate(y=adopted_sarc_alt2)

#tests for the within group comparisons, Only MiM
tests_wncs_mim <- data %>%
  select(cs,pd,y) %>%
  filter(cs=='MiM') %>% 
  group_by(cs) %>% 
  t_test(y ~ pd, detailed = TRUE)
#gotta add the y position...
tests_wncs_mim <- tests_wncs_mim %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values and estimates nicely 
tests_wncs_mim <- tests_wncs_mim %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs_mim)

#tests for the within group comparisons, excluding MiM
tests_wncs <- data %>%
  select(cs,pd,y) %>%
  filter(cs!='MiM') %>% 
  group_by(cs) %>% 
  t_test(y ~ pd, detailed = TRUE)
#gotta add the y position...
tests_wncs <- tests_wncs %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values nicely 
tests_wncs <- tests_wncs %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs)

#tests for the between group comparisons, but excluding MiM
tests_bwcs <- data %>%
  select(cs,pd,y) %>%
  filter(cs!='MiM') %>% 
  group_by(pd) %>%
  t_test(y ~ cs, detailed = TRUE)
#gotta add the y position...
tests_bwcs <- tests_bwcs %>% 
  add_xy_position(x = "cs" , group = "pd" , dodge = 0.8) #fun = "mean_sd" ,
#format the p values nicely 
tests_bwcs <- tests_bwcs %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_bwcs)

#tests for the diffs-in-diffs
tests_did_lm <- data %>% 
  select(cs,pd,y) %>%
  filter(cs!='MiM') %>% 
  lm_robust(y ~ pd * cs , . , se_type = "stata")
print(summary(tests_did_lm))
didb <- summary(tests_did_lm)$coefficients["pdDiscounted:csSDM", "Estimate"] %>% 
  round(., digits=3)
didp <- summary(tests_did_lm)$coefficients["pdDiscounted:csSDM", "Pr(>|t|)"] %>% 
  round(., digits=3)
#might have to make a separate object with x y and the values we want
did <- data.frame(a = paste("DiD = ",didb," ","p = ",didp))
print(did)



ggbarplot(data, x = "cs" , y = "y", 
          fill = "pd", 
          add = c("mean_ci"),
          position = position_dodge(0.8),
          xlab = "Counselling style" , 
          ylab = "Fraction of clients",
          title = "Fraction of clients who adopted a SARC - Broad definition", 
          ylim = c(0,1),
          label = TRUE,
          lab.nb.digits = 3,
          lab.hjust = -.35,
          lab.size = 3) +
  #add in the p-values and differences, all default set at 1 so have to manually lower them down
  stat_pvalue_manual(tests_wncs, label = "d = {est}; p = {prounded}",
                     tip.length = 0.01,
                     bracket.nudge.y = -0.35,
                     label.size = 3) +
  stat_pvalue_manual(tests_wncs_mim, label = "d = {est}; p = {prounded}",
                     tip.length = 0.01,
                     bracket.nudge.y = -0.35,
                     label.size = 3) +
  stat_pvalue_manual(tests_bwcs, 
                     label = "d = {est}; p = {prounded}",
                     color = "black", 
                     tip.length = 0.01, 
                     step.increase = 0.08,
                     bracket.nudge.y = -0.225,
                     label.size = 3) +
  #add the DID estimates 
  geom_text(data = did, 
            aes(x = 2.5,  y = .975, label = a),
            size = 3) +
  annotate("segment", x = 1.8, xend = 3.2, y = .95, yend = .95, size=0.25) +
  annotate("segment", x = 1.8, xend = 1.8, y = .95, yend = .94, size=0.25) +
  annotate("segment", x = 3.2, xend = 3.2, y = .95, yend = .94, size=0.25) +
  #make it look nice
  theme_bw() +
  scale_fill_manual(values = met.brewer("Egypt", 2)) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12),
        legend.text = element_text(size=12),
        legend.position= "bottom",
        legend.title = element_text(size=12),
        panel.grid.minor = element_blank(), 
        panel.grid.major = element_blank(),
        plot.caption = element_text(hjust=0,size=10)) +
  guides(fill=guide_legend(title="LARC prices:")) 
#  labs(caption = "SARCs include: Pill (COC or POP), Injectable, LAM, Condoms, Calendar+Cycle Methods, Abstinence, Emerg. contr.")

ggsave("output/figures/FigS3-view-prices-sarcs-broad.png", dpi = 320, scale = 0.8, height = 16, width = 24, units = "cm")
